Celestin Apprentice 7/Apprentice-Release7.iso/Source Code/Pascal/Applications/NIH Image 1.62b11/src/LeastSquares.p

Text File | 1994-11-30 | 7KB | 319 lines

unit LeastSquares;
{Contains the curve fitting routines based on the Simplex method described in the article " Fitting Curves to Data "}
{in the May 1984 issue of Byte magazine, pages 340-362.}

interface

    uses
        QuickDraw, Palettes, Printing, globals, Utilities, graphics;


    const
        nvpp = 2;
        maxn = 6;
        maxnp = 20;
        alpha = 1.0;
        beta = 0.5;
        gamma = 2.0;
        lw = 5;
        root2 = 1.414214;
        MaxError = 1e-7;

    type
        ColumnVector = array[1..maxnp] of extended;

        vector = array[1..maxn] of extended;
        datarow = array[1..nvpp] of extended;
        index = 0..255;


    var
        m, n: integer;
        done: boolean;
        maxx, maxy: extended;
        i, j: index;
        h, l: array[1..maxn] of index;
        np, npmax, niter, maxiter: integer;
        next, center, smean, error, maxerr, p, q, step: vector;
        simp: array[1..maxn] of vector;
        data: array[1..maxnp] of datarow;
        filename, newname: string;
        yoffset: integer;


    procedure DoSimplexFit (nStandards, nCoefficients: integer; xdata, ydata: ColumnVector; var Coefficients: CoefficientArray; var residuals: ColumnVector);


implementation


    function f (p: vector; d: datarow): extended;
        var
            x, y, ex: extended;
    begin
        x := d[1];
        case info^.fit of
            StraightLine:
                f := p[1] + p[2] * x;
            Poly2:
                f := p[1] + p[2] * x + p[3] * x * x;
            Poly3:
                f := p[1] + p[2] * x + p[3] * x * x + p[4] * x * x * x;
            Poly4:
                f := p[1] + p[2] * x + p[3] * x * x + p[4] * x * x * x + p[5] * x * x * x * x;
            ExpoFit:
                f := p[1] * exp(p[2] * x);
            PowerFit:
                if x = 0.0 then
                    f := 0.0
                else
                    f := p[1] * exp(p[2] * ln(x)); {y=ax^b}
            LogFit:
                begin
                    if x = 0.0 then
                        x := 0.5;
                    f := p[1] * ln(p[2] * x)
                end;
            RodbardFit:
                begin
                    if x = 0.0 then
                        ex := 0.0
                    else
                        ex := exp(ln(x / p[3]) * p[2]);
                    y := p[1] - p[4];
                    y := y / (1 + ex);
                    f := y + p[4];
                end; {Rodbard fit}
        end; {case}
    end;


    procedure order;
        var
            i, j: index;
    begin
        for j := 1 to n do
            begin
                for i := 1 to n do
                    begin
                        if simp[i, j] < simp[l[j], j] then
                            l[j] := i;
                        if simp[i, j] > simp[h[j], j] then
                            h[j] := i
                    end
            end
    end;


    procedure sum_of_residuals (var x: vector);

        var
            i: index;
    begin
        x[n] := 0.0;
        for i := 1 to np do
            x[n] := x[n] + sqr(f(x, data[i]) - data[i, 2])
    end;


    procedure Initialize;
        var
            i, j: index;
            firstx, firsty, lastx, lasty, xmean, ymean, slope, yintercept: extended;
    begin
        firstx := data[1, 1];
        firsty := data[1, 2];
        lastx := data[np, 1];
        lasty := data[np, 2];
        xmean := (firstx + lastx) / 2.0;
        ymean := (firsty + lasty) / 2.0;
        if (lastx - firstx) <> 0.0 then
            slope := (lasty - firsty) / (lastx - firstx)
        else
            slope := 1.0;
        yintercept := firsty - slope * firstx;
        case info^.fit of
            StraightLine:
                begin
                    simp[1, 1] := yintercept;
                    simp[1, 2] := slope;
                end;
            Poly2:
                begin
                    simp[1, 1] := yintercept;
                    simp[1, 2] := slope;
                    simp[1, 3] := 0.0;
                end;
            Poly3:
                begin
                    simp[1, 1] := yintercept;
                    simp[1, 2] := slope;
                    simp[1, 3] := 0.0;
                    simp[1, 4] := 0.0;
                end;
            Poly4:
                begin
                    simp[1, 1] := yintercept;
                    simp[1, 2] := slope;
                    simp[1, 3] := 0.0;
                    simp[1, 4] := 0.0;
                    simp[1, 5] := 0.0;
                end;
            ExpoFit:
                begin
                    simp[1, 1] := 0.1;
                    simp[1, 2] := 0.01;
                end;
            PowerFit:
                begin
                    simp[1, 1] := 0.0;
                    simp[1, 2] := 1.0;
                end;
            LogFit:
                begin
                    simp[1, 1] := 0.5;
                    simp[1, 2] := 0.05;
                end;
            RodbardFit:
                begin
                    simp[1, 1] := firsty;
                    simp[1, 2] := 1.0;
                    simp[1, 3] := xmean;
                    simp[1, 4] := lasty;
                end;
        end;
        maxiter := 100 * m * m;
        n := m + 1;
        for i := 1 to m do
            begin
                step[i] := simp[1, i] / 2.0;
                if step[i] = 0.0 then
                    step[i] := 0.01;
            end;
        for i := 1 to n do
            maxerr[i] := MaxError;
        sum_of_residuals(simp[1]);
        for i := 1 to m do
            begin
                p[i] := step[i] * (sqrt(n) + m - 1) / (m * root2);
                q[i] := step[i] * (sqrt(n) - 1) / (m * root2)
            end;
        for i := 2 to n do
            begin
                for j := 1 to m do
                    simp[i, j] := simp[i - 1, j] + q[j];
                simp[i, i - 1] := simp[i, i - 1] + p[i - 1];
                sum_of_residuals(simp[i])
            end;
        for i := 1 to n do
            begin
                l[i] := 1;
                h[i] := 1
            end;
        order;
        maxx := 255;
    end;


    procedure new_vertex;
        var
            i: index;
    begin
        for i := 1 to n do
            simp[h[n], i] := next[i]
    end;


    procedure DoSimplexFit (nStandards, nCoefficients: integer; xdata, ydata: ColumnVector; var Coefficients: CoefficientArray; var residuals: ColumnVector);
        var
            i, j: integer;
            d: datarow;
    begin
        np := nStandards;
        m := nCoefficients;
        for i := 1 to np do
            begin
                data[i, 1] := xdata[i];
                data[i, 2] := ydata[i];
            end;
        Initialize;
        niter := 0;
        repeat
            done := true;
            niter := succ(niter);
            for i := 1 to n do
                center[i] := 0.0;
            for i := 1 to n do
                if i <> h[n] then
                    for j := 1 to m do
                        center[j] := center[j] + simp[i, j];
            for i := 1 to n do
                begin
                    center[i] := center[i] / m;
                    next[i] := (1.0 + alpha) * center[i] - alpha * simp[h[n], i]
                end;
            sum_of_residuals(next);
            if next[n] <= simp[l[n], n] then
                begin
                    new_vertex;
                    for i := 1 to m do
                        next[i] := gamma * simp[h[n], i] + (1.0 - gamma) * center[i];
                    sum_of_residuals(next);
                    if next[n] <= simp[l[n], n] then
                        new_vertex
                end
            else
                begin
                    if next[n] <= simp[h[n], n] then
                        new_vertex
                    else
                        begin
                            for i := 1 to m do
                                next[i] := beta * simp[h[n], i] + (1.0 - beta) * center[i];
                            sum_of_residuals(next);
                            if (next[n] <= simp[h[n], n]) then
                                new_vertex
                            else
                                begin
                                    for i := 1 to n do
                                        begin
                                            for j := 1 to m do
                                                simp[i, j] := (simp[i, j] + simp[l[n], j]) * beta;
                                            sum_of_residuals(simp[i])
                                        end
                                end
                        end
                end;
            order;
            for j := 1 to n do
                begin
                    if (simp[h[j], j] - simp[l[j], j]) = 0 then
                        error[j] := 0
                    else if simp[h[j], j] <> 0 then
                        error[j] := (simp[h[j], j] - simp[l[j], j]) / simp[h[j], j]
                    else
                        error[j] := (simp[h[j], j] - simp[l[j], j]) / simp[l[j], j];
                    if done then
                        if abs(error[j]) > maxerr[j] then
                            done := false
                end;
        until (done or (niter = maxiter));
        ShowMessage(concat('interations=', long2str(niter), crStr, 'max interations=', long2str(maxiter)));
        for i := 1 to n do
            begin
                smean[i] := 0;
                for j := 1 to n do
                    smean[i] := smean[i] + simp[j, i];
                smean[i] := smean[i] / n;
            end;
        for i := 1 to m do
            Coefficients[i] := smean[i];
        for i := 1 to nstandards do
            begin
                d[1] := xdata[i];
                Residuals[i] := ydata[i] - f(smean, d);
            end;
    end;


end.